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Abstract 

During the coalescence of binary black holes, gravitational waves carry linear momentum away 
QQ . from the source, which results in the recoil of the center of mass. Using the Effective One Body 

approach, that includes nonperturbative resummed estimates for the damping and conservative 

(N 

' parts of the compact binary dynamics, we compute the recoil during the late inspiral and the 

!>■ : 

subsequent plunge of non-spinning black holes of comparable masses moving in quasi-circular orbits. 
Further, using a prescription that smoothly connects the plunge phase to a perturbed single black 
hole, we obtain an estimate for the total recoil associated with the binary black hole coalescence. 
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^ • We show that the crucial physical feature which determines the magnitude of the terminal recoil 

^H^. is the presence of a "burst" of linear momentum flux emitted slightly before coalescence. When 

using the most natural expression for the linear momentum flux during the plunge, together with 
^ ' a Taylor-expanded (f/c)^ correction factor, we find that the maximum value of the terminal recoil 

is ~ 74 km/s and occurs for r] = ^^]^_^yi — 0.2, i.e., for a mass ratio m2/mi ~ 0.38. Away from 
this optimal mass ratio, the recoil velocity decreases approximately proportionally to the scaling 
function /(r/) = r/^ (1 — 4r/)^^^ (1.0912 — l.OAr] + 2.92?]^). We comment, however, on the fact that 
the above 'best bet estimate' is subject to strong uncertainties because the location and amplitude 
of the crucial peak of linear momentum flux happens at a moment during the plunge where most of 
the simplifying analytical assumptions underlying the Effective One Body approach are no longer 
justified. Changing the analytical way of estimating the linear momentum flux, we find maximum 
recoils that range between 49 and 172 km/s. 

PACS numbers: 04.30Db, 04.25.Nx, 04.80.Nn, 95.55.Ym 
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I. INTRODUCTION 



During the coalescence of compact binaries, along with energy and angular momentum, 
the system radiates linear momentum. The loss of linear momentum via gravitational radi- 
ation results in the recoil of the center of mass of the binary. It is astrophysically important 
and desirable to obtain a dependable estimate for the velocity of the center of mass of 
comparable- mass black hole binaries undergoing coalescence 0,Q]. Notably, in models of 
massive black hole formation involving successive mergers, recoils large enough to eject co- 
alescing black holes from dwarf galaxies or globular clusters would effectively terminate the 
process. This motivation has recently led several authors to estimate the recoil velocity of 
coalescing black hole binaries by means of various methods QflQ. Ref. 3 employed black 
hole perturbation theory to describe the motion of a test mass moving in a black hole back- 
ground, i.e., the case where the symmetric mass ratio t] = („™'j|_^^)a 1- They combined a 
numerical estimate of the recoil velocity up to the Last Stable Orbit (LSO), with two crude 
estimates for the recoil acquired during the subsequent plunge phase. Then they assumed 
that their test-mass estimates could be proportionally scaled up to the comparable-mass 
cases (4?7 ~ 1) with the function 

/(r/) = r/2^/r34^, (1) 

which appears as an overall factor in the leading ("Newtonian") analytical estimate of 
the recoil, as first computed in Ref. 0. The final estimates of Ref. P] range (for non- 
spinning black holes) between a lower value ~ 54 {f {f]) / fmax) km/s and an upper one 
~ 465 (/(r/)//max) km/s, where = /(r/^ax) = /(0.2) = 0.0178885. Another set of 

estimates was obtained from an approach that employs a mixture of numerical relativ- 
ity and blacky hole perturbation theory for the merger of comparable mass non-spinning 
In particular, for a mass ratio m2/mi = 0.5, corresponding to a value 



black holes 



Q. 



rj = 0.2222, close to the value ?7inax = 0.2, where /(?]), given by Eq. (^, reaches its maxi- 
mum value [/(0.2222)//(0.2) = 0.01646/0.01789 ~ 0.92], Ref. ^ estimates a recoil velocity 
~ 250 ± 150 km/s. Finally, using an analytical estimate, which is further discussed below, a 
maximum recoil (reached for r] ~ 0.2 ) equal to 250 ± 50 km/s was obtained in Ref. 

Summarizing, the recent estimates are consistent with a maximum recoil velocity ~ 
250 km/s for non-spinning black holes In contrast, we shall estimate here, by 

using the Effective One Body (EOB) approach to binary black hole dynamics, detailed in 



Refs. 
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ll|, that the maximum recoil velocity for non-spinning coalescing black 
holes is probably significantly smaller, and of the order ~ 50 — 74km/s. However, we shall 
conclude that this estimate is rather uncertain because it depends on a specific way of 
describing the linear momentum fiux during a crucial phase of the plunge which is (mildly) 
relativistic, and has not been yet analytically studied in detail. When changing our preferred 
asssumptions for describing the linear momentum fiux, we find maximal recoil velocities that 
vary in the range 49 — 172 km/s. 

Let us recall that the coalescence of isolated black hole binaries may be viewed as con- 
sisting of the following three phases. The first phase is that of gravitational-radiation-driven 
slow inspiral in quasi-circular orbits. This leads, after a very long period of adiabatic shrink- 
age of the orbital separation, to the binary approaching its LSO. At the LSO, the inspiral 
phase changes to some sort of plunge, dominated by general relativistic strong-field effects. 
This plunge phase results in the merger of the two black holes and the dynamics of the final 
black hole, thus formed, may be described in terms of black hole quasi-normal modes. 

During the early inspiral, the dynamics of compact binaries can be described very accu- 
rately by the post-Newtonian (PN) approximation to general relativity. The PN approx- 
imation to general relativity allows one to express the equations of motion for a compact 
binary as corrections to Newtonian equations of motion in powers of {v/cY ~ GM / {(?R), 
where f , M, and R are the characteristic orbital velocity, the total mass and the typical 
orbital separation of the binary respectively. For a compact binary, treated to consist of 
non-spinning point masses, of arbitrary mass ratio rj = where mi and m2 are the 

masses of the components and M = mi + m2, the leading "Newtonian" contribution ^ to 
linear momentum fiux JFI and the associated instantaneous recoil v* i.e., the instanta- 



neous velocity of center of mass can be derived from the investigations that dea^ 
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t either with 



15|. This 



wave-zone fiux computations, or with near- zone radiation-reaction ones 
allowed Ref. to derive the following explicit leading order results (for a circularized orbit 
of radius R) 



^ Or rather "quasi-Newtonian", as this leading contribution corresponds to 3.5PN ©(w^/c'')] radiation 
reaction terms in the equations of motion for the binary. 
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where j(r]) denotes the combination introduced in Eq. (^. The function ]{r]) vanishes 
both in the hmit of extreme mass ratio (test-mass hmit, 77 ^ 0), and in the case of equal 
mass binaries (77 = 1/4). It reaches the maximum value /max = 0.0178885 for ?7max = 0.2, 
corresponding to a mass ratio ma/mi = (3 ± Vs) /2 ~ 0.3820 or 2.618. In the following, 
we shall generally (though not systematically) employ units where G = c = \. 

The above leading-order results do not allow one to reliably estimate the final recoil 
of a coalescing binary black hole. Indeed, on the one hand they neglect higher-order PN 
corrections that might become fractionally large near the LSO and during the subsequent 
plunge, and on the other hand, they do not take into account the crucial transition between 
,ua.-*cu.a. .sp.a. and plunge. A. exe.pUfied bv p.eviou. „o..s fl B Q, one can 
think of several different ways of going beyond the above results, given by Eqs. As 
first attempted by Detweiler and Fitchett one can use perturbation theory around 
black hole backgrounds to calculate the linear momentum flux emitted by a test particle 
moving around a black hole, assuming then a scaling proportional to jjj]) in order to go 
from the test-mass case (r/ ^ 1) to the comparable- mass one. Ref. [la| considered only 
circular motions (both above and below LSO), while Ref. combined information about 
circular orbits above the LSO with crude estimates of the linear momentum flux during 
the plunge following the LSO crossing. However, we note that this approach is yet to be 
applied for the relevant case of (non-geodesic) plunging motions. The Lazarus program, 
which employs a mixture of perturbation theory and numerical relativity, to study the late 
inspiral and the merger of black hole binaries is not limited to the case 77 <^ 1. However, 

n 

these simulations are limited to rather short evolution time spans. For instance, Ref. |^ 
mentions that the simulations are accurate only for "less than 15 M in time", which is 
significantly smaller than the 3PN-accurate estimate of the orbital period of a comparable- 
mass black hole binary at the LSO: — ^ P]- This forces state-of-the-art numerical 
simulations to employ initial data sets where the two holes are already quite near to each 
other (and, in fact, probably closer than the LSO). It is indeed a challenge to construct initial 
data sets for tight black hole binaries, and which do correspond to the physically correct 



"no-incoming radiation" criteria. [See Refs. 17,^^ for prescriptions to construct 'realistic 
initial data sets'.] Presence of spurious radiation in the initial data sets is likely to dominate 
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and thereby invalidate the estimate of the hnear momentum flux. However, a comparison 
between numerically constructed initial data sets and analytically determined characteristics 
of tight binary systems, performed in Ref. have indicated the ability of the Effective 
One Body (EOB) approach to describe rather accurately the initial data obtained 

by the helical Killing vector method 2^ , which is carefully aimed at minimizing the amount 
of spurious radiation present in the initial state. 

In view of the above situation, we consider that our present "best bet" for obtaining a 
dependable estimate for the gravitational recoil during the late inspiral and plunge phases of 
a black hole binary consists in employing the EOB approach. This approach has already been 
used to compute complete gravitational waveforms emitted during the inspiral and merger of 
(non-spinning and spinning) black holes 8, 11|. As first proposed in Ref. ^, this is achieved 
by considering the waveform emitted by the binary beyond the LSO, through the subsequent 
"plunge", down to, approximately, the "light ring" {R ~ 3M), and by matching it there to 
a "ring down" signal constructed using the quasi-normal modes of the resultant final black 
hole. As this (extended) EOB approach will be central to the present paper, let us recall 
the main arguments of Ref. ^ for proposing the apparently bold strategy of analytically 
describing the plunge beyond the LSO, down to i? ~ 3 M. A first argument is that the EOB 
approach is a resummation technique which was carefully devised to work not with the 
usually considered, badly convergent, PN-expanded equations of motion or flux quantities, 
but instead with a (EOB) re-summed Hamiltonian and a (Fade) re-summed damping force 
showing^o sign of bad behavior during most of the "plunge" . In particular, it was found 
in Ref. (sf that the word "plunge" to qualify the dynamics beyond the LSO is a misnomer, 
and that this phase is better thought of as being still a quasi-circular inspiral motion, even 
down to the light ring i? ^ 3 M. Indeed, it was found that the quasi-circularity condition 
{R <^ R(f) remains satisfied with good accuracy beyond the LSO, down to i? ~ 3 M. This 
is illustrated in Fig. Q] which shows (for rj = 0.2) the evolution of the "azimuthal" {g'^'^p^) 
and "radial" ( g^^pji) kinetic energies during the plunge down to the light-ring, i? ~ 3 M. 
The crucial point is that the ratio TZ = g^^p%/g''^'^p^ stays significantly smaller than one 
during the entire plunge. Its value at i? ~ 3 M is T^-ught-ring — 0.281. 

As for the idea of matching the gravitational wave emission to a quasi-normal-mode 
(QNM) "ring-down" signal around i? ~ 3 M, let us recall that it was realized long ago that 
the basic physical reason underlying the presence of a QNM-type merger signal, that ends 
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the plunge signal, was that the / > 2 gravitational waves emitted by the collapsing system 
are strongly filtered by the potential barrier, centered around i? ^ 3 M, describing the radial 



propagation of the gravitational waves 
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Recently, Blanchet, Qusailah and Will [5| have employed an "approximation" to the 
EOB method in the sense that they 'assume that the plunge can be viewed as that of a "test 
particle" of mass fi = moving in the fixed Schwarzschild geometry of a body of mass 

M'. They also assumed that the effect of radiation reaction damping on the plunge orbit 
may be ignored. They then matched, in various ways, a circular orbit at the Schwarzschild 
LSO, i.e., 6 M, to a suitable plunge orbit. By contrast, we shall use here the full EOB 
approach 0|, which does not need to "assume" that the plunge can be viewed as that of 
suitable test particle, but instead proves it [see Refs. 3, Ol], and which does not need 
to match a circular orbit to a plunge orbit at the LSO, because it automatically embodies 
a smooth transition between the "inspiral" and the "plunge". Let us also note that the 
"effective test body" used in the EOB method does not evolve in a fixed Schwarzschild 
geometry of mass M, but instead in a deformed Schwarzschild background, whose geometry 
was algorithmically derived to 2PN accuracy in Ref. [3], and to 3PN order in Ref. [see 
also Ref. for the incorporation of spin effects]. In addition, while Ref. 0] formally let 
their "test particle" fall down to the horizon at i? = 2 M, an important ingredient of our 
approach will be to match the plunge signal to a QNM-based ring-down one at i? ~ 3 M. 

On the other hand, an important result of Ref. concerns the higher-order PN correc- 
tions to the "Newtonian" linear momentum flux, given by Eq. Using the multi-polar 



)ost-Minkowskian approach 



3, 31, 



271 . l28l . 129^ , and its higher-order implementations 



Ref. 



has gone beyond the previous IPN-accurate studies of recoil effects, 
available in Ref. |33||, by including both the 1.5PN order "tail" contribution and the next 

n 

2PN order corrections. Ref. p finds that the linear momentum flux at infinity, for binary 
systems in circular orbits, is given by a PN expansion of the form 

464 



j(BQW) 



105 



(3) 



^ For the tcst-particlc case, this foUows from the expUcit form of the Rcgge- Wheeler- Zcrilh effective poten- 
tials. In the comparable-mass case, we must contemplate the I > 2 binary gravitational wave signal as 
propagating in the spacetime generated by the binary system, and approximate the latter (when the two 
holes are closer than 3 M, and when considering the waves propagating in the domain 7? > 3 M) by the 
external geometry of a single hole (of mass-energy ~ A/). 
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where /(r;) is given by Eq. above ^, and 

v^^{Mnf\ (4) 

with Vt = dif/dT denoting the orbital angular velocity. The factor F{vi^; 77) yields the 2PN- 
accurate "Taylor-expanded" PN-corrections to the linear momentum flux (when the latter 
is expressed in terms of the above defined v^) and it reads 

FK; r/) = 1 + F^iv) vl + vl + F,{r^) , (5) 

with 

^ , , 452 1139 

= -— -^r/, (6a) 

FM = ^vr, (6b) 

71345 36761 147101 , 

^4(77) = \ 71 H 71^ . (6c) 

^" 22968 2088 ' 68904 ' ^ ^ 

Finally, A* in Eq. Q is a tangential unit vector directed in the same sense as the relative 
orbital velocity v"^ = v\ — v^, v\ and being the velocities of the masses mi and m2 
respectively. Note that the test-mass limit (77 0) of the function F{v^) has been first 
numerically evaluated in Ref. Q]. As we shall see below, one of the important differences 
between our treatment and the one of Ref. will concern the continuation of the linear 
momentum flux Eq. Q (derived for circular orbits above the LSO) to the (non circular) 
plunging orbit below the LSO. 

In the next section, we present our prescription to compute the linear momentum flux 
and the related velocity of center of mass. Section ITTTl contains a summary of the 'modified' 
EOB approach that is used to describe the late stages of binary inspiral and plunge, followed 
by a detailed account of the numerical procedure that will result in the determination of the 
associated gravitational radiation driven recoil. We also present in that section analytical 
insights into our numerical estimates. In Sec. IIV[ we describe how we smoothly match the 
merger and the resultant ring down phases and the computation of the recoil of the final 
black hole. We present our results, conclusions and future directions in Sec. El 



^ We conventionally assume henceforth that mi > m2 so that = +^1 — 4i]. 
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II. QUASI-NEWTONIAN FORMULAS FOR LINEAR MOMENTUM FLUX AND 
RELATED RECOIL 



In the EOB formalism, one finds that the relative orbital dynamics of a binary black hole 
system is most conveniently described in a "Schwarzschild-like" coordinate system, to which 
is associated an "effective metric" of the form 

dsls = -A{R) dT^s + dR^ + (de^ + sin^ 6 d^^) . (7) 

We shall work here to 2PN accuracy, in which case the "effective metric coefficients" A{R) 
and D{R) are given by 

2M /mV 



DiR) = l-6r/ (^-J . (8b) 

It was shown in Ref. that the complicated and badly convergent second post-Newtonian 
expanded dynamics of a binary system could be mapped onto the much simpler (and better 
convergent) dynamics of an auxiliary test particle falling along a geodesic of the effective 
metric, Eq. (jZj). Note that, even in the equal mass limit (77 = |), the effective metric 
coefficients, given by Eqs. (jH}, differ only slightly from those of a Schwarzschild metric (i.e 
As{R) = 1 — ^^;Ds{R) = 1). As emphasized in Ref. 0| this property makes it useful 
to describe the EOB dynamics in the Schwarzschild- type coordinates of Eq. ((7j), rather 
than, say, in Arnowitt-Deser-Misner or harmonic coordinates which would lead either to a 
(badly convergent) infinite series of PN corrections, or to more complicated "resummed" 
expressions. 

One of the important features of the present study will be to express the flux of linear mo- 
mentum radiated away from a compact binary directly in terms of the quasi- Schwarzschild 
coordinates R and (f used in the effective one body metric, Eq. ((Tj), and of their time- 
derivatives, notably the angular velocity Q = Our work will often rely on the use of quan- 
tities having simple (and "quasi-Newtonian") expressions in terms of quasi- Schwarzschild 
coordinates R and (p. Let us first motivate this use of quasi-Newtonian quantities expressed 
in quasi-Schwarzschild coordinates. 

In the test-mass limit, it is a striking feature of Schwarzschild coordinates that they of- 
ten allow one to convert Newtonian results into exact, or near-exact, Einsteinian results. 
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A famous example of that is the location ( a la Mitchell-Laplace) of the Schwarzschild 
horizon which is correctly given by using the purely Newtonian energy conservation: 
c(r)^/2 — GM/R = c(oo)^. In addition, the angular frequency along circular geodesies 
in a Schwarzschild background is described, in Schwarzschild coordinates, by the usual 
Kepler law: GM = Q^R^, so that the linear velocity v = QR is given by the usual New- 
tonian formula = GM/R. Here we shall use the remarkable fact that this closeness 
extends to gravitational radiation properties. In particular, the total energy flux emitted 
by circular geodesies into gravitational waves is numerically very well approximated by 
the simple quasi-Newtonian formula obtained by writing the leading-order quadrupole for- 
mula [dE/dt oc {d^Iij/ dt^Y oc f2^(Jjj)^] in Schwarzschild coordinates. Indeed, this yields 
dE/dt = ^rfVL^I^ which is quite close to the complete general relativistic answer 2^: even 
at the LSO, -Rlso = 6M, the quasi-Newtonian result ^rfVL^R^ is only 12% smaller than the 
full Einsteinian one, and the agreement is better for orbits above the LSO. A look at Fig. 2 in 
Ref. shows that a similar type of agreement holds also for the flux of linear momentum 
down to the LSO. Note that is crucial in this comparison to "interpret" quasi-Newtonian 
results in terms of Schwarzschild coordinates. For instance, if one were to insert harmonic 
coordinates Rh in the quadrupolar result dE/dt = ^tfVL^R^ (and uses a corresponding 
harmonic-coordinate Kepler law M = ^"^Rf/) one would obtain an estimate for the energy 
flux which would be larger than the correct one, at the LSO {Rh = 5), by a factor ~ 2.19. 

Regarding the plunging orbits below the LSO, another important feature of our treatment 
is that we do not wish to insert in the leading-order, "quasi-Newtonian" , energy and linear 
momentum fluxes the usually assumed "Kepler-type" law relating the angular velocity Q 
to the radius R. Indeed, Kepler's law (which reads GM = Q^R^ in the test-mass limit 
?7 ^ and in Schwarzschild coordinates) is only valid, below the LSO, along the physically 
irrelevant sequence of unstable circular orbits corresponding to a maximum of the effective 
radial potential^. The "violation" of Kepler's law during the plunge will be illustrated in 
Fig. 2 below. 

With this motivation, let us derive from scratch the quasi-Newtonian result for the linear 
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Indeed, this maximum of the effective potential corresponds to an unphysical angular momentum > 
plf''^. In contrast, the physically relevant plunge motion corresponds to < p]f''-' and thereby to a 
particle gliding down a flatfish effective potential having no maximum (nor minimum) anymore, i.e, a 
potential near but below the lowest radial potential plotted in Fig. 1 of Ref. 
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momentum flux. We start from tlie following leading order formula, available in Ref. [1 

1 r (4) (3) (4) (3) (3) >, 

336 TT 1^ J 
where J-^ and jFp are the x and y components of the linear momentum flux. Here /'"^ and 
gim fignote the "mass" and "spin" (or "current" ) radiative multipole moments of the binary, 

(3) 

while /'™ denotes the third time derivative of /'™. Under complex conjugation P"^ and 
transform as 

We display below the relevant P"^ and S^"^ required to compute the leading order contri- 
bution to jFp + jFp for compact binaries in circular orbits, taken from Ref. j35[ , 
7-22 2 921 o 

— = -VlO^rjR'e-''^, _ = __v/l0^r/v/^^/?3^]e-^^ (11a) 

j31 2 /33 2 

where = The exact expression for, say, P^(x — ^ — ^ contains several terms 
proportional to ^ In the following, we consider an inspiralling and plunging 

relative orbit. For such an orbit, the derivatives do not vanish. However, as already 
mentioned above, it was pointed out in Ref. that even during the "plunge" following 
the LSO crossing, the radial motion, characterized by remained small compared to 
the azimuthal one R We shall take advantage of this fact to simplify the expression 
of the time-differentiated multipole moments entering Eq. Q by keeping only the terms 
proportional to the time-derivatives of the azimuthal angle ip. We neglect also = ^ 
compared to ^ = fi. This yields the simplified expression 

( + I order = -l^ f{n) ^'e'^. (12) 



An important difference between expression ((T2j) and the earlier quoted expressions for linear 
momentum fluxes, namely, Eqs. Q and (jH)), is that the proportionality to R^ VP was directly 
obtained from the original flux formula, Eq. ©, without explicitly using any "Kepler-like 
equation" linking R to VL. We shall see later that this difference significantly affects the 
estimate of the final recoil velocity associated to the linear momentum flux (fT^ . 

In order to obtain the velocity of the center of mass, we then invoke linear momentum 
balance, namely, 

M ^ «om + ^ <m) = - (-^^ + ^ ^^) , (13) 
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where v^^^ and v"^^^ are the x and y components of the center of mass velocity vector Vcom- 
At this stage, it is convenient to introduce rescaled, dimensionless radial, time and fre- 
quency variables, namely, r = -0, t = and lo = ^ = MVl. This leads to the following 
differential equation for v^^^ + i v^^^ 

^ «om + ^ vy^J = z — f{v)r'u;'e^^. (14) 

This leading order "quasi-Newtonian" result will be the basis of our investigation. We 
shall also discuss below how to use the 2PN correction terms derived in Ref. [see Eq. Q] 
to improve the accuracy of Eq. (fT^ . 

In this paper, in order to obtain an estimate for the velocity of the center of mass during 
the late inspiral and subsequent plunge phases, we shall numerically integrate Eq. (fT^ along 
with the differential equations that define the EOB dynamics. 

In the next section, we summarize the EOB dynamics applicable to non-spinning compact 
binaries of arbitrary mass ratio moving in quasi-circular orbits during the inspiral phase . 
We also describe, in some detail, how we solve the relevant set of differential equations to ob- 
tain an EOB based estimate for the recoil during the late inspiral and the subsequent plunge 
phases. We shall also complement our numerical estimates by analytic arguments allowing 
one to understand in simple terms the main characteristics, and the order of magnitude, of 
our results. 



III. THE LATE INSPIRAL, PLUNGE PHASES AND THE ASSOCIATED RE- 
COIL USING THE EOB APPROACH 

Let us first summarize the EOB approach relevant for describing the inspiral and plunge 
phases of a compact binary. At the 2PN accuracy, the mapping between the full two-body 
2PN dynamics, and the much simpler geodesic dynamics in the EOB metric, given by Eq. ((7j), 
leads to an EOB dynamics described by the following Hamiltonian (expressed in the scaled 
variables r = jj,t=^,uj = ^ = M^l, and in polar coordinates) 

n{r,pr,p^) = - 
V 

where A{r) and B{r) [see Eqs. ((Tj), and (jH)) above] are given by 

A(r) = l-^ + ^, (16a) 



l + 2r] 



'A{r) 1 + 



Pr 



B(r] 



(15) 



11 



B{r) ^£!fl='(l-^4]. (16b) 



A{r) A{r) \ 

More precisely, the explicit form of the EOB equations of motion read 

dr _ dH{r, p-,.,Py>) 
dt dpr ' 



(17a) 



dt dp^ 
dpr dn{r,pr,p^) 

(17c) 



dt dr 

. (17d) 



dp^ 



dt 

The right-hand side of the last equation expresses the loss of angular momentum under 
gravitational radiation reaction. Its explicit form will be discussed below. 

As mentioned earlier, to obtain an estimate for the recoil during the late inspiral and 
subsequent plunge, we solve along with the above set of differential equations, the one for 
^com + ^^com' given by Eq. dUD, namely, 

d 464 

^ «om + ^ ^cV) = ^ /(^) - F . (18) 

Here the supplementary factor F = 1 + OixP'^ + Oiv^^ + 0(t'^), resulting from Ref. j^, 
is added to improve the accuracy of the leading-order, quasi-Newtonian result, given by 
Eq. to the 2PN level. Its explicit form along our quasi-circular, sub LSO, orbits is 
discussed in the following subsection. 



A. Inclusion of 2PN corrections in the fluxes of linear and angular momenta 

In this subsection, we describe the construction of 2PN accurate expressions for the recoil 
(linear momentum flux) factor F in Eq. (fH^ . as well as for the radiation reaction force Tip 
(the angular momentum flux), appearing in Eq. ()17d|) . 

Let us start by discussing the value of the correcting factor F (and of its analog in the 
energy flux) during the adiabatic inspiral phase. During this phase, our construction is 
facilitated by the fact that the orbital dynamics closely follows the one parameter sequence 
of stable circular orbits that exists above the LSO. In the EOB formalism, these orbits 
represent the minima, with respect to r, of the Hamiltonian HdTcir,p^) = H{r,pr = 0,p^). 
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Equivalently, we see from Eq. ()15p that they are obtained by minimizing with respect to r 
the effective potential 

w{r,p^) = A{r) + , (19) 

with A{r) given by Eq. ()16aj) . Minimizing w{r,Pip) with respect to r yields the following 
relation linking r to [see Eqs. (4.5) and (4.6) of Ref. ^] 

P;icirc = r ^_^~^^^ (20) 

Inserting the latter result in the definition of the angular velocity, namely oo = 

also considered along circular orbits [i.e., Uchc = ^"^"g^^'^'^'' ], then yields a relation connecting 

cUcirc to r. This 2PN generalization of Kepler's third law reads [see Eqs. (4.8) of Ref. j^] 

'''^- = Kl + 2!y(v^-l))- ^^^^ 
In the test mass limit, ?7 — > 0, we recover the well-known fact that circular orbits in a 
Schwarzschild geometry (in Schwarzschild coordinates) satisfy the standard Kepler law: 
(^2^3 _ jg f^i^QT^ traditional to use as PN order parameter = cu^^'^ = 0{v/c), 

or equivalently x^^ = v"^ = iJ^I"^ = 0{v'^/c^), to describe all possible PN corrections, be 
they proportional to the square of the linear azimuthal velocity v^p = ur, or to the grav- 
itational potential m = Indeed, when 77 = 0, we have the simple, Kepler-like links: 

/ \2 1 2 

and to the sub-LSO quasi-circular orbits ), let us introduce the function 



— ^ — "'^ — ' To extend these simple links to the comparable mass case rj 



^(r. P^) = ^ 3;^ , (22) 



l + 2r][ ^Jw{r,p^) - 1 

1 _ 3i? 

and the definition 

r^=r\^{r,pp)\ . (23) 

These definitions are such that, along circular orbits, we can still write a simple Kepler- 
looking law 

^V3 = l, (24) 

as well as its usual consequences, such as [ur^)^ ~ V~ ~ ^"^^^ = "^S- then use these 

relations to rewrite any 2PN-accurate result expressed (along circular orbits) in terms of 
Voj = co^^^ in terms of u, r and if). 
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For instance, the 2PN accurate linear momentum flux p, Eq. Q, is proportional to 
v^j^ F{v^) = uj'^ r^jF^ur^) = lJ^ {r ^tl)^!'^^ Fiior 't\}^l'^\ Our approach leads us to considering 
that the basic "quasi-Newtonian" expression for the linear momentum flux is proportional to 
^7j,5 i^ggg j^q_ (ji2j) above]. In other words, we are naturally led to writing the 2PN-accurate 
flux in the form of Eq. (|18p with a 2 PN- correct ion factor F given by 

F{r,p^) = U(r,p^) j F{urij'/') . (25) 

Let us flrst note that for circular orbits above the LSO (for which all the above reasonings 
are fully justifled) the "correcting factors" linked to the function tp are very close to 1. More 
precisely, if we consider the case t] = 0.2 (which is the most important one) tends to 
1 when r oo, and as r decreases ip flrst decreases to reach a minimum ipmm — 0.9882 
around r ~ 9.2. Afterwards, it increases to reach ■j/'lso — 0.9921 when r = tlso — 5.8. 

1 /3 

Note that the factor V^lso — 0.9974 modifying the azimuthal velocity v^^, = lu r in Eq. ()25|) 
differs only by ~ 3 x 10"^ from unity. As for the total 2PN correcting factor F, one can see 
that it represents, above the LSO, a relatively modest modiflcation of the quasi-Newtonian 
momentum flux. If we evaluate it by inserting the straightforward 2PN-expanded version of 
the function F, Eq. into Eq. ()25|1 . we flnd a result of order 1.24 at the LSO. 

Note also that, in the EOB approach, it is natural to consider as basic PN-ordering 
parameter the azimuthal velocity, 

= u;r , (26) 

which is an invariantly defined quantity ^. 

Up to this stage we have been assuming that we were considering quasi-circular orbits 
corresponding to a local minimum of the effective radial potential. This happens when one 
is above the LSO. In contrast, when considering the continuation of the orbit below the 
LSO, the circular orbits [and their consequences, such as Eq. (PT|) ] are no longer physically 
relevant because they correspond to unstable maxima of w{r,p^), given by Eq. (fT^ . When 
considering quasi-circular orbits below the LSO, one should, in principle, re-derive from 
scratch the 2PN-accurate linear momentum flux, without assuming any Kepler-like law of 
type Eq. (PT|) or Eq. (j2H). The 2PN corrections then become functions of three independent 

^ Indeed, in the EOB formalism, both u) and the effective-metric Schwarzschild radius r are invariant 
quantities. 
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variables: r, p<^, and Pr, or equivalently, r, u, and pr. As pointed out in Refs. 0,0] the motion 
remains "quasi-circular" during the plunge in the sense that the contributions linked to p^ 
stay numerically small compared to those linked to p^, allowing one to neglect p^. The 2PN 
corrections during the plunge then become functions of two independent variables, namely, r 
and Pip. In view of the arguments recalled above pointing to a remarkable closeness between 
exact Einsteinian results and quasi-Newtonian results expressed in terms of Schwarzschild- 
type coordinates, we consider it likely that the momentum flux during the plunge follows 
more or less the quasi- Newtonian behavior JF oc r^cu^. To ensure continuity with the 2PN- 
correcting factor Eq. (j2Sl), which is present above the LSO, we shall assume here that 2PN 
corrections below the LSO are sufficiently well estimated by continuing to use the expression 

Note that this assumption differs from the one made in Refs. and P| which consisted 
in continuing to use the expressions giving the 2PN corrections as functions of = cu^^^. 
Within the spirit of the EOB formalism, we feel that it is not very plausible to continue 
to use as basic PN ordering parameter, and to express quantities only in terms of it. 
Indeed, the definition of v^^ makes sense only so far as a Kepler-like law relating u to r 
continues to hold. This is no longer the case below the LSO (as will be illustrated in Fig. 
2 below). In absence of such a Kepler law, we prefer to remain close to what is suggested 
by the leading-order quasi-Newtonian result JF oc r^cj^. We shall further discuss below the 
importance of this choice. 

Let us apply the same philosophy to the estimate of the radiation reaction term JF^,, 
appearing in Eq. (jl7d|) . below the LSO. It was shown in Refs. and [S| that a good 
estimate for JF^ above the LSO is given by 



^-= = -|,.„:Miki5), (27) 



"^pole 



where /dis is a 2.5PN (f ^-accurate) Pade approximant for the angular momentum flux. It 
is deflned, e.g, in Eqs. (3.28)-(3.36) of Ref. [Note, however that the 2.5PN coefficient 



^ Note, however, that we no longer assume the Unk, provided by Eq. (|20|l . Indeed, this hnk exhibits an 
infinite growth of as r tends to the (jy-modified ) "hght-ring" , where r'^ — 3 7'^ + 5 ij — 0. As wc know 
instead that stays below its LSO value during the plunge (and evolves much more slowly than r), it 
seems a priori better to express the PN corrections only in terms of well behaved quantities, such as r, w, 
and p^. 
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there must be corrected to a new value due to Ref. [37[.]. The definition of /dis depends 
on the choice for the location of the "pole" fpoie- Following Ref. js^ [and also Ref. jsl], we 
use the value of fpoie defined in Eq. (3.37) of Ref. j^. [As shown in Ref. 2^ the precise 
choice of Wpoipis not very important and, for instance, the value Vpoie = "^^poic — 73 "would 
suffice.] Ref. j8| proposed to continue using Eq. (jTTj) . expressed in terms of v^^ = cu^^^, even 
below the LSO. Here, consistently with the arguments presented above, we shall instead 
use a different continuation for JF^ below the LSO. To derive it we need to know what is 
the analog, for the angular momentum loss, for the "quasi- Newtonian" result, displayed in 
Eq. 0. Consistently with what was briefly mentioned above about the energy flux, we 



know [see, e.g., Eq. (4.23) of Ref. [15|] that the leading term in the angular momentum 

(2)_ (3) 

loss ^ is oc X /^^. Remembering the leading-order expression for the quadrupole 
moment J^^, given in Eqs. we see that the "quasi-Newtonian" expression for JF^ i 



IS 



OC — ^jj^ — X — dr3 therefore (in the quasi-circular approximation r <^ r (p) 

oc r^u;^. This shows that one should rewrite the leading factor in Eq. ()27p as i.e., 
= uj^ r^. In other words, this leads us to using, below the LSO, the following expression 
for the radiation reaction force JF^ in Eq. (j2Z|), 

"^polc 

where the factor is a. function of r and p^, deflned in Eq. ()22|) and where u = uj{r,pr,Pip) 
is deflned by the second equation in Eqs. (fTTjl . 

Let us flnally discuss the question of the re-summation of the 2PN-accurate correction 
factor, Eq. ()25|1 . When comparing the straightforward, PN-expanded versions of the 2PN 
factor F{v^), given by Eq. (0) and derived in Ref. j^, entering the 2PN accurate linear 
momentum flux to its analog in the energy (or angular momentum) flux, namely, 

Fe{v^) = Taylor = 1 + F^.iv) vl + F^,{r^) vl + < + (29) 

^polc 

where FE2 = -(^ + ir/),FE3 = 47r, and Fe4 = ("w + w ^ + § [see Ref. Q], 
one notices that these two "Taylor expansions", i.e., expressions for F{v^^) and F-£,{Vi^), are 
rather similar. The corresponding Taylor coefficients F„ [in ^^d [in JFe or Tip] 

have the same signs, similar sensitivities to the value of 77, and roughly similar magnitudes. 
Indeed, we can roughly consider that F„ ~ 1.3 FEn- In addition, the same argument which 
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was used in Ref. j36[ to show that the "exact" function Fe{v) ( analytically continued from its 
behavior above the LSO) has, in the limit r/ — > 0, a pole at Wp~j|! = ( i.e. at the light ring ) 
can be applied to the function F{v), appearing in Eq. (jS)) for jFp^^^-*, to conclude that F{v) 
also has, when 77 ^ 0, a pole at the same (light-ring value) v2~^ = This suggests that a 
Pade re-summation of the F{v), given by Eq. (0), might improve the convergence behavior of 
F{v), which is currently known only up to 2PN accuracy. On the other hand, the similarity 
between the two Taylor expansions -F„ and F^n (i-e. the fact that F„ ~ 1.3 Feu) suggests that 
both the successive Taylor and Fade approximants of F{v) will have convergence properties 
similar to the corresponding approximants of Fe- By looking at the convergence properties^ 
(when ?7 — s> 0) of the Taylor or Pade approximants of Fg, as displayed in Fig. 3 of Ref. js^, 
we observe that, among the approximants of FN order v"' with n < 5, the best one is the 
f ^-accurate (2.5FN level) Fade approximant [i.e the one used in Eq. (jTTj) above]. However, 
the 2.5FN coefficient, F^lrj), is not currently known for the analogous linear-momentum flux 

^*(BQW) p^^^ p.g 3 j^gf ^ 

we expect the f ^-accurate Taylor approximant of F{v) 
to overestimate F'^^'^'^^{v) when v < 0.4 and to underestimate it when v > 0.4 (note that 
vlf'~'{ri = 0.2) ^ 0.414). On the other hand, we expect the f ^-accurate Fade approximant to 
F{v) to follow F'^^'^^{v) better, but to underestimate it for all values of v. In the absence of 
knowledge about the 2.5FN contributions -^5(77), we shall compare here the results obtained 
both from using the f ^-accurate Taylor approximant, given by Eq. (0), and its corresponding 
Fade approximant, of the form 

FA.) = ^, (30) 
where we take, for simplicity, fpoic = and with 

9iv) = , , \. . (31) 



^ Note that the argument used in Ref. |3 (namely: the closeness of the v^- and w^-accurate Taylor approx- 
imants suggests that a good convergence is reached with 2PN accuracy) is not conclusive in view of what 
happens for the similar Taylor expansion of Fe{v). Indeed, one can easily check (say, when 77 = or 
77 = 0.2) that the v'^ and v'^ accurate Taylor approximants of Fe{v) are close to each other, while the next 
w^-accurate, Taylor approximant of Fe{v) is quite far away from both of them (and also from the exact 
result, when ?/ ^ 0). This is linked to the fact, emphasized in Ref. [s^ . that Taylor approximants have 
rather erratic convergence properties as the PN order increases (while, by contrast, the Pade approximants 
have a more monotonic convergence, though they tend to accumulate somewhat below the exact result). 
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Note that g{v) is constructed quite similarly to fDis{v), i.e. by applying 
Eqs. 13.29), (3.31), (3.34) and (3.35) of Ref. Q, while replacing Eqs. (3.32) and (3.33) in 
Ref. js| by the coefficients F2, F^, and F4 entering Eq. ^ above. 



B. Initial conditions for the dynamics 

Before we present our numerical results, let us explain how we prescribe initial conditions 
when solving these differential equations. The initial conditions for r and if are arbitrary 
and prescribe the initial radial separation of the binary in the center of mass frame and its 
associated angular position. The initial values for pr and are obtained with the help of 
the adiabatic approximation to the EOB inspiral, introduced in Sec. IV (A) of Ref. j^. This 
approximation is obtained by imposing pr = in the EOB dynamics, which implies that the 
effective body follows an adiabatic sequence of circular orbits with decreasing energy due to 
the emission of gravitational radiation. This zeroth-order adiabatic approximation (which 
turns out to be enough for our purpose) provides the following expressions for p^ and u 

.2 \ 1/2 



^ I adiab 



r — 3 77) 
r-^ — 3 + 5 

(1 - ^ 



r3 [1 + 27] ^{ 



(32b) 



where z{r) = with A{r) given by Eq. (jl6a|l . 

The initial value for v^^^ + i v^^^^ is obtained by using 

^com linitial = ^ V 1 " 4 7/ , <m |initial = ^ V 1 " 4 , (33) 

derivable from Eqs. (0). In the next subsection, we present analytical insights into the 
physical behavior underlying our numerical estimates. 



C. Linear momentum loss during inspiral and plunge 

Before numerically implementing our strategy for estimating the recoil during the late 
inspiral, the subsequent plunge and the final merger, let us outline the main physical features 
of our calculation*. Let us first note that the final recoil velocity is essentially given by an 

® See below for the effects of ring-down. 
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integral of the form 

/ + 00 
dta(t)e^-W, (34) 
■oo 

where (p{t) is the orbital phase, while the "amplitude" a{t) = |jFp| is proportional to r^u'^ F. 
The first important point to realize is that the value of the above integral X is dominated 
by what happens in the time intervals where the amplitude a{t) varies in a "non-adiabatic" 
manner. Indeed, let us first assume that the amplitude a{t) always varies in an adiabatic 
manner with respect the orbital phase (p{t), i.e let ^ -C 0. This can be conveniently for- 
malized by replacing the phase factor e**^ by e'^'^^'^, where e is a formal "small parameter" 
measuring how small the ratio ^/^ ~ ^{^) is- Using the fact that ip{t) is a monotonic 
function of time, one can use ip, instead of t, in the above integral. After replacing a{t) by 
A{ip) = g|, we get X(e) = J^^ dipA{ip) e'^^'^^ . Using e'^/^ = f (^^J^) and integrating 
by parts, we find X(e) = ie J d(pA'{(p) e^'^^^, where we employed the vanishing of A{(p) at 
±00 (see below). Repeating this procedure, we get 

/+00 
rfy^^(")(^)e'^/% (35) 
■00 

where A^"'\'p) is the n-th derivative of A w.r.t (p. This (well-known) result means that, 
when e ^ 0, X(e) = (9(e") for any integer n. In other words, X(e) vanishes faster than any 
power of e, if a{t) varies adiabatically during the whole process. In most cases, this means 
that X(e) is exponentially small ~ e"'"/'^^, and therefore numerically negligible compared to 
the naive estimates of the type X ~ ctmax At, where At is a characteristic variation time, say 
At ~ ^ with LJ = (p, that one might have been tempted to make. 

This mathematical reminder shows that the actual magnitude of the momentum flux 
|jFp| = a(t) during the inspiral and the plunge is secondary with respect to the question of 
knowing the characteristic time-scale on which |jFp| varies during the plunge. If ^ stayed 
always small compared to the orbital frequency u = (f, the recoil would be a non-pertubative 
effect; and it would be practically hopeless to try to estimate it by starting from approximate 
analytical expressions. [On the other hand, we would know that the recoil is exponentially 
small, so that it would be astrophysically negligible.] However, the study of the time evo- 
lution of the amplitude a{t) during the EOB plunge shows that, while it remains adiabatic 

<^ (f) during most of the inspiral and plunge, it becomes barely non-adiabatic near the 
moment where a{t) reaches a maximum [at which point, the criterion for non-adiabaticity 
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must involve the second time derivative of a(t)]. This is illustrated in Fig. |2l which shows 
the evolution of the magnitude of the momentum flux, |^p|, during the late inspiral and 
the plunge. This figure makes it clear that the characteristic evolution time scale for a{t) 
is shortest near its maximum, i.e. for a (scaled) radius Tmax — 3.5. Before discussing the 
consequences of this fact, let us outline how one can analytically understand why |jFp| has 
the behavior exhibited in Fig. |21 

The main factor determining the behavior of |jFp| is the product cj^ in Eq. p4|) . During 
the plunge, the evolution of r and uj are governed by the EOB equations of motion, namely 
Eqs. (ITTj) . In these equations, the radiative damping "force" JF^ is crucial to drive the slow 
inspiral and to trigger the plunge, but was found to play a minor role once the plunge is well 
on its way As a consequence, the evolution of uj during the plunge is approximately 

given hj 00 = with ^ constant as well as 7i ^ constant ( zero-damping approx- 

imation). In this approximation, one thereby finds that uo is approximately proportional to 
the ratio so that we can write, during the plunge, the approximate link 

U;plunge(r)^..LSO^^j^4^. (36) 

The change in behavior of uj{r) between the inspiral (where Kepler's third law uo'^r^ ^ 
constant holds approximately), and the plunge [where Eq. (jHU)) holds instead] is illustrated 
in Fig. El 

As a consequence of Eq. (j36|) . we find that the behavior of the linear momentum flux 
during the plunge is approximately given by 

|^p|cxr=.'<xi(l-i + ?|)', (37) 



Denoting m = ^, it is easily seen that |jFp| oc (1 — 2 m + 2rju^y reaches a maximum value 
when (9 — 32 u^ax + 60 r/uj^^x) — 0- When rj = 0.2 ( which corresponds to the maximum of 
the overall factor /(?]), given by Eq. (P), and thereby approximately to the maximum recoil 
), this gives Umax — 0.29044, corresponding to ^'^ax'^*"^ = — 3.4431. This analytical 
argument agrees well with our numerical results. Indeed, we find that \J-'p\, computed using 
the quasi-Newtonian version of Eq. ()18|1 . has a maximum around r^^'^'^^'^^^ — 3.501. One can 
even go further and analytically study the behavior of |JFp| near its maximum. The most 
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important time scale there is defined by the curvature of \J^-p\ near its maximum: 

^max 

''max = / w ^ \ max ; (^^) 



where, for notational simphcity, we henceforth denote by jFp the modulus of the momentum 
flux. 

An important dimensionless quantity associated to the time scale Tmax is the "quality 
factor" Q = coimaxTmax associated to the "resonance peak" of JFp. Indeed, values of Q of 
order unity mean (as we shall find) that the evolution of jFp near its maximum is just fast 
enough to be non-adiabatic there. 

In view of the discussion above, this means that the recoil, i.e. (modulo a factor i) the 

integral X, given by Eq. (j34p . will be dominated by what happens near the maximum of 

jFp. Therefore, we can analytically estimate X by replacing a(t) = J^p{= |^p|) by the local 

approximation ( with i =t — t^^^ ) 

1 ffl'^T \ '^^^ / 1 f2 \ f2 

a{t) ^ + - ( P = JT'^ h - - -^J ^ JT'^e 

and the phase ip{t) by 



(39) 



- 1 - 1 

'^{t) — V^max + 0max ^ + 2 '^^^^ ^ ^ ^ 2 ^ ' ^^^^ 

If we then consider the recoil acquired up to some given time ?, its is given by a truncated 
Gaussian integral 

Km + ^^Jin^^-^r'^e^^- r rfte-^"*"+^*, (41) 

^—00 

where (posing e^ax = c^max r^ax ) 

1 1 

a = — i Cjmax = (1 - i Cmax) , /? = « ^^max ■ (42) 

^max ^max 

When i' gets positive and large with respect to Tmax ( so that i' is effectively ~ +00 ), we 
can estimate the total integrated recoil using the standard complex Gaussian integral: 

»+oo 1^ 

dye-^'^y'+^'y = V2^ . (43) 



This yields, for the modulus of the corresponding total integrated recoil 

^max 
(1 + 4ax) 



2 2 

T 1 '-^max ^max 
' ma">f — 77 



^com ' '' ^coml — v^'i^p ^ xiM \ / 
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This approximate analytical result vividly illustrates the preceding discussion. Indeed, if the 
evolution of T-pit) were adiabatic {Q = tUmaxTmax ^ 1 ) the total integrated recoil would be 
exponentially small ( even if jFp^^ gets large ). 

We have already indicated above, see Eq. (j37j) . how one can analytically determine 
the location on the r-axis of the maximum of jFp. By analytically expanding Eq. (|37p 
around its maximum, one can also get an analytical expression for the product Tmax^^max- 
As the reasoning above also gave the variation with r of the angular frequency, namely, 
Lj{r)/uLso ~ [-^{^)/''^'^]/[-^{''^LSo)/''^Lso]' '^^^ obtain analytical estimates of uj-^ax and of 
^max/''^max- Finally, to get analytical estimates of all the quantities entering Eq. (jUj) we 
need an analytical estimate of r around r = rmax- This can be obtained (though only with 
modest accuracy) by using again the zero-damping approximation to write that the energy 
is approximately conserved during the plunge: hence 



A{r) 



l+pl/B{r)+pl/r' 



wlso , (45) 



with Wlso ~ ^ij'Lso)i)- + P^lso / ''^lso) ■ -^l- approximately determines the value of pr 
during the plunge. From it, one then deduces the value of r by using the Hamilton equation, 
Eq. ()17ap . See Fig. 4 of Ref. ^ for a plot of r during the plunge. 

These analytical approximations allows one to obtain estimates for all the quantities en- 
tering the crucial Eq. (j44j) . and thereby to obtain an analytical estimate of the expected total 
recoil velocity Vcom- We found that the results agrees within a few percent with the numbers 
one can extract from our full numerical simulations. The complete set of relevant quantities, 
extracted from our simulations, for the behavior around the time where |jFp| reaches its max- 
imum value are (for r] = 0.2): r^ax — 3.501, fmax — —0.113, ciJmax — 0.1255, cj^ax — 2.952 x 
10-^ iJ^lmax ^ 2.039 X 10-\ |J^p|inax ^ 3.683 X 10"^ r ~ 7.440, Q ~ 0.9337, e - 0.1634, and 
^ 1.007. Of particular importance is the value of the quality factor, namely, Q ^ 0.9337. 
The fact that it is of order unity means that a net integrated recoil is acquired soon after 
T-p{t) reaches its maximum value, i.e. soon after the plunge has fallen below r ~ 3.501, and 
therefore before reaching the light ring radius rir = 3. 

Finally, we can analytically estimate, using Eq. ()44|). the final recoil that one might 
expect. We find 

Kom + ^ ^^c^omP"^^^^^^^' ^ 74.06 F..ax km/s , (46) 

where we recall that /(0.2) = 0.0178885, the maximum value reached by /(r/) (when r) = 
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0. 2), and where -Fmax denotes the value of the 2PN correction factor at r = rmax — 3.501. 
By definition, the quasi- Newtonian estimate corresponds to taking F = 1. 

IV. TRANSITION FROM PLUNGE TO RING-DOWN, AND GRAVITATIONAL 
RECOIL DURING RING-DOWN 

Though the analytical estimate, given by Eq. (j33I), is interesting by the physical informa- 
tion it conveys [effect dominated by the maximum of J-p{t), dependence on f{ri) and Fmax, 
and the obtainment of a small pure number from high powers of numbers "of order unity"], 
let us hasten to add that it is only an approximation to the real terminal recoil. Indeed, the 
above estimate, was obtained by taking the formal limit t +oo in the truncated Gaussian 
integral, Eq. ()41|) . However, as we have already indicated in the introduction, the physics 
behind the approximate analytical formulas, Eqs. (fT^. ((HI), or (fTHj) . changes when r reaches 
the "light ring" r ^ 3. Following the analogous estimate of complete waveforms in Ref. j^, 
we propose here to estimate the contribution to the recoil due to the merger of black holes 
by formally terminating the plunge when the scaled radial coordinate gets around r ~ 3, and 
by matching there the relevant time derivatives of the radiative multipole moments during 
the late plunge phase to corresponding "ring-down" multipole moments, constructed from 
appropriate quasi-normal mode contributions. 

Let us first discuss why it is important to ensure as smooth a matching as possible during 
the transition from plunge to ring-down. To see this, let us consider again the approximate 
form of the final recoil velocity, given by Eq. ()34|1 . but let us now divide the full time interval 
in two phases: an inspiral -|- plunge phase, lasting from t = — oo up to some tmatch, followed 
by a ring-down phase, lasting from tmatch up to t = +oo. The total recoil will be the sum of 
two contributions of the form ^ 

dt (O'(t) G ) Iplungc ) -^ving / dt ('^(t) G ^ Iring • l^"^) 

We now focus on the contribution to Xtot = Xpiungc + Xnng that is formally linked to any 
"mismatch" between the two behaviors of the linear momentum flux around t = tmatchj 

1. e. to any discontinuity between i ^j^'^p^^^'^s.c _|_ ^ j^wiungcj _ j-^^^^-j 6* ]piunge' considered for 

^ Actually, the integral during the ring-down phase is a sum of terms oc e+ "^ and e" "''. See discussion 
below. 
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t < Wtch, and i ( jri^"'^s 



[a{t) considered for t > tmatch- The effect of 



p -r t ^ p 

any discontinuity around t = tmatch can be obtained by summing the "edge contributions" of 
the two semi-infinite integrals, Eq. ()47|). i.e. the contributions hnked to the upper or lower 
cut-off t = tmatch- Thesc "cdgc contributions" have been worked out in Ref. to next- 
to-leading order in "adiabaticity expansion" [i.e. in powers of the formal small parameter 
e, introduced in Eq. (j34|) above, by replacing e**^ by e**^/^ ], by using the "integration by 
parts" technique introduced above for showing that, in the absence of any discontinuity, the 
integral X(e) = J dta{t) e*'^^*)/^ vanishes faster than any power of e. Adding two terms of 



the type of Eq. (3.17) in Ref. 



38[ , the total edge contribution is of the form 



plunge ring 



a{t) 



iip{t) 



plunge 



ring 



(48) 



where the square bracket 



- plunge 



^(a(t),(^(t),0(t). 



ring 



on the right-hand side of the 

above equation denotes the difference, J^(apiungc(tmatch), V^piunge(Wtcii), 0piungo(imatch), •••)- 
J^(aring(Wtcii),V5ring(Wtcii),0ring(Wtcii), •••)• This analytical resuh highlights the following 
fact: any discontinuity between the amplitude, the phase, or any of their time-derivatives 
across tmatch "will contribute to the final recoil velocity. Therefore, if we want to minimize 
the spurious effects linked to our describing the smooth transition between the plunge and 
the merger by a fictitious sharp transition happening at t = tmatch, we should try to match 
as many derivatives as possible of jFp oc a{t) e**^*^*^ across t = tmatch- On the other hand, we 
are going to see that, even after having matched as well as possible JFp across t = tmatcin 
there remains a (non-spurious) "edge" contribution linked to the physical change of behavior 

across t = tmatch- 

To see this, let us consider in more detail how one can implement a physically motivated 
matching across t = tmatch ( corresponding to rmatch — 3). All the physical effects which are 
important for the present study ( flux of energy related to J^^, and flux of linear momentum. 



jFp) can be expressed as integrals over a sphere at infinity with integrands proportional to 
the local gravitational wave energy flux = TjS^ oc ( r hij ) , where hij is the TT-gauge 



dimensionless gravitational wave amplitude, and hij is its time derivative [see, for e.g., 
Ref. Qj]. This motivates us to try to match as well as possible the quantity r hij{t, r, 6, 0), 
where 6 and are polar angles on the sphere at infinity, between plunge and ring-down. The 
"radiative multipole moments" that enter the multipole expansion of r hij are, by definition, 
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the (/ + l)-th time derivatives of the l-th mass {P"^) and spin (or current) (S*'™") multipole 

moments It is therefore most natural to match P"^ and 5*'™ across t = tmatch- For the 
evaluation of jFp at the leading order, the relevant radiative moments, as seen in Eq. Q, 

(3) (3) (4) (4) 

are P^^; S"^^^; I^^^ and I^^^. These terms correspond to gravitational waves, emitted by 
the two black holes, of multipolarity: ( / = 2, m = ±2, even parity); ( / = 2, m = ±1, odd 
parity); ( / = 3, m = ±1, even parity); ( / = 3, m = ±3, even parity), respectively. As a first 
approximation we can consider that these gravitational waves propagate (for radii larger 
than the radial distance r{t) separating the two black holes) on a Schwarzschild background, 
of mass Ms = -^[0^' ~ M = mi + m2, approximately representing the (physical) spacetime 
outside the two holes. Therefore, when r{t) gets smaller than about 3, the relevant modes 
fij^p will be strongly filtered by the corresponding Regge-Wheeler-Zerilli effective potential 
Vil^'^^^"^'^\r). This filtering can be approximated by saying that, when the source of a mode 
(/,m, tt), where vr denotes the parity, falls below r = 3, the corresponding outgoing wave 
mode can be described by a superposition of quasi-normal modes (QNM's) of the same 
multipolarity (/,m, vr). 

Several nice simplifying features of gravitational wave propagation on a Schwarzschild 
background are that: (i) the effective potential Vi^{r) does not depend on the "magnetic 
quantum number", m, (ii) Vi"{r) is real, and (iii) though VJ'^™'^ ^ V","'^'^, they have the same 



spectrum of QNM complex frequencies [for a review of QNM's see Ref. |4l|]. For each value 
of the multipolar order /, there is a double infinite sequence of QNM complex frequencies, 

say 

a^^ = ain±i ujin , (49) 

where n = 0, 1, 2, ...,and and uiin are both real and positive [so that a'^ = (o"^^)*]- The 
notation here is that the n-th QNM mode belonging to the multipolarity / decays, when 
t — > +00, proportionally to e"'^'"* = e"""'"* e'''*'^'"*. For each value of / the fundamental 
QNM mode = is the least-damped one, i.e. the one with the smallest value for 
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For simplicity, we use here the nomenclature of Ref. . In the multipolar post-Minkowskian formalism, 

Refs. |2iEiEiE3,EiEi0ElEi,the "radiative" moments are defined as J7'™ -/'™ and y'™ -5'", 
i.e. as the moments entering the multipole expansion of rhij. In the latter nomenclature, the moments 
that most directly enter the quantities that we need would be C/'™ and (and would include all required 
'tail' effects). 

We leave to future work the refinement consisting in using modes propagating over a Kerr background. 
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Finally, our matching procedure consists in joining, as smoothly as possible, across t = 

(3) (3) (4) 

tmatch, each relevant multipolar mode entering rh^^p, namely, ll^l^^it) , Sl^^^^^{t) , I^^l^^{t) , 

(4) 

and -^plunge (^)' obtained for t < tmatch by differentiating Eqs. (fTTj) in the quasi-circular ap- 
proximation (r ^ r(j)), to corresponding "ring down" multipole moments, made of sum 
of decaying QNM modes. For instance, this leads (after scaling out the total mass M) to 
matching 

(3) 16 

5 



plunge it)=^ — VT0^vr{tY u{tf e-2*^(*) (for t < Wh) , (50) 



where r{t),ip{t),uj{t) = ip{t) are obtained by numerically integrating the EOB dynamics, 
Eqs. f)17|) . to a corresponding "ring down" radiative moment of the form 

^nig (t) = E {C'n (^'') ^"'^^"^ + ^nin e"'^-^} (for T^t- > 0) , (51) 

n=0,l,.. 

where o"^, n = 0,1,.., are the QNM frequencies, Eq. (jlH)), belonging to the multipolarity 
/ = 2, and where Cj^(/^^) denotes, for each n, two independent complex coefficients. Indeed, 

(3) 

/^^ (t) being complex, there are no reality conditions relating C^(/^^) and C~(/^^) ( in spite 
of the fact that o"^ and o"^ are related by complex conjugation). 

If we include in Eq. ()51|) only the first two complex conjugated fundamental QNM modes, 

(3) 

crj) and (J20, we observe that 1^^^ (t) contains two arbitrary complex coefficients Cq^P"^) 
and Cq{P^). These two complex coefficients can be chosen so as to ensure not only that 

(3) (3) 

^plunge (t = ^match) agrccs with I^^^ {t = t^^tch) = (I^^) + Cq (P^) , but also that the 

(3) (3) 

(numerically computed^^) time derivative — s^^^e — agrees, when t = tmatch with — — = 
-^2+ C+{P') - a^-Q C^iP'). This yields 
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/ON (3) 
^- 7-22 (J.\ _|_ -'plunge 

"20 -'plunge y'^J ^ dt 



Coil 



22\ 



^20 - ^: 

(3) 

"^20 -^plunge (^) + dt 



(3) ^j22 
22 ( + \ \_ plunge'- 



C^20 



t^t^^ , (52a) 



For a smooth match, one should no longer use the quasi-circular approximation, r <C r (/j, when computing 

(3) 

the time derivatives of /piu,-,g(, {t). 
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(3) (4) 

Similarly, we can match, in a once-differentiable (C^) manner, Spl^^^^^t) , r^l^^^^{t) , and 

(4) 

-^plunge (^) ^'^ ^i'^S down moments of the form 

(3) 

^ nng it) = C^iS'') e-^tor + (S'') c'^^- ^ (53a) 

(4) 

/ring (^) = C+il'') e--^-- + e-'^ao- , (53b) 

(4) 

4fng (^) = C^{I''') e"'^^"^ + Co-(/'') e-'^ao- . (53c) 

Each pair of complex coefficients C^{A4) is then given as a linear combination of 
A4piunge(^match) and plunge (^match) of the type, given by Eq. §^ above. Finally, we 
can use the complex conjugation relations, Eq. to match, in a manner, the remain- 

(3) (4) 

ing required radiative moments P~'^{t) and J^~^(t) entering Eq. This does not introduce 
new, independent coefficients as 

i-rC^il'"') = [{I'"")*] = [Co^(/'™)]* . (54) 

Note also that to match the multipole moments entering the leading order linear momentum 
flux, we need to know only two conjugate pairs of complex QNM frequencies, namely from 
Refs. and H, 

afo = 0.08896 ± 0.37367 i, (55a) 
a± = 0.09270 ± 0.59944 i. (55b) 

Having so determined continuations of the various relevant multipole moments during the 
merger phase, we get an estimate of the flnal recoil, in the leading-order (quasi-Newtonian) 
approximation by integrating, from — oo to tmatch and then from tmatch to +oo, the linear 
momentum balance equation 

^ 1 r (4) (3) (4) (3) (3) 



dt 



^con. + ^<J = -^|v^/^-^/^^ J^^/^-^ -14z P-^S^'Y (56) 



Here the "radiative moments " /'"^ and S*'™ appearing on the right-hand side (RHS) are 

(i+i) 

given by: (i) when t < thatch by Eq. and similar "plunge moments" /^igc "^pSnge 

obtained by differentiating ( while neglecting r <C r 0) Eqs. (Ilip . and (ii) when t > tmatch by 

(i+i) (i+i) 

analytical QNM-based "ring-down moments" ^ring, >S''™g, deflned by Eqs. above. The 
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continuity of the moments entering the RHS of Eq. ()56|) ensures that the hnear momentum 
flux jFp [defined by the RHS of Eq. is continuous, as well as its first derivative, across 

t ^match' 

As explained above this matching ensures that one did not introduce leading-order 
spurious contributions linked to edge effects. At the same time, this matching procedure 
generically introduces discontinuities in the second time derivative of jFp. We then see from 
Eq. ()48|) that there will be sub-leading spurious contributions linked to such discontinuities 



in To study the eventual numerical importance of these higher-order edge effects, we 

have also implemented an improved matching procedure consisting of including, for each 
radiative multipole moment, the first two conjugate pairs of QNMs in Eq. (|5ip. i.e. both 
n = ( fundamental QNMs) and n = 1 (first excited QNMs). The new required QNM 



frequencies, available in Refs. 
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42], are 



af^ = 0.27391 ±i 0.34671, (57a) 
af^ = 0.28130 ±i 0.58264. (57b) 

As the QNM sums, Eq. ()51|1 . now include 4 arbitrary complex coefficients Cq ,Cq ,Ci ,Ci , 

(3) 

we can uniquely determine them by demanding that each radiative moments ( say ), 

,(3) 

together with their first three numerically computed time derivatives ( for j < 3) match 

across tmatch- 

The matching procedure presented so far was based on considering the leading-order, 
quasi- Newtonian, expression, Eq. 0, for the fiux of linear momentum. When considering 
the 2PN correction factor F, as in Eq. (jl8|) . we should, in principle, both include more 
multipolarities in Eq. Q, and PN-corrections in the expressions for individual radiative 
moments, Eqs. (fTT|) or Eq. (j^Uj) . As we found (see below) that contributions to the recoil due 
to the ring-down phase are relatively small, we decided, for simplicity, to use a less-rigorous, 
but much simpler, 2PN-level matching procedure. The procedure we used consisted in 
continuing to use the leading-order fiux, as in Eq. (jH)), but to "improve" the "brick radiative 

(3) (3) (4) 

moments", J^"*, S*^™, P"^, it contains by multiplying each of them by a factor v F, e.g. we 
modify Eq. (|Hnjl to 

/„\ improved 

(3) 16 



^''plunge =^VF-x/T0^r/r(t)2a;(t)3e-2^^W (for t < Wh) • (58) 
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Then we match each of these "improved plunge moments" to a corresponding "ring-down" 
one, given by a QNM sum of the form, given by Eqs. (jM]) . Again this matching can be done 
in a Ci (2 QNM's) or ^ 4 QNM's) manner. 

V. RESULTS 

Having presented our methodology, let us now discuss the results that we obtained, and 
their interpretation. Let us consider first the leading-order, quasi-Newtonian approximation, 
i.e. Eqs. (jl2|) and (jl4j) . together with the leading-order (2 QNM's per moment) matching 
to the ring-down phase. We plot in Fig. ^ for the case of 77 = 0.2, the magnitude of 
the linear momentum flux \Tp\, together with their two separate components jFp and Tp, 
as functions of time. The maximum of |JFp| is reached for t = t^^^ ^ 4149.50 (which 
corresponds to Tmax — 3.501, in an evolution for which the initial separation when t = 
was r = 15, while the matching to the ring-down phase was done at t = tmatch — 4153.50, 
which corresponds to rmatch — 2.933. Note the rather fast (and oscillatory) decay of the 
individual components of jFp during the ring-down. Indeed, we see from Eqs. ()5H) to ()53|) 
that, during the ring-down, jFp is a sum of contributions proportional either to e~v'^2o+<^3oJ = 
g- («2o +030 ) T g-j{e2 1^20 +1:3 1^30 )t (-where €2 = 1 = 63) or to e~v °"2o j ^ _ ^_2Q,20Tg-j(e2+e2)'^20T 
(where = 1 = e'l). From Eqs. and (fS^. we see that the slowest exponential decay 
is oc e~^"2'''^, which decays on a characteristic time scale Tring = ~ 5.62. Though this 
is significantly smaller than the orbital period near the LSO, note, however, that this time 
scale is comparable both to the characteristic time scale for the variation of J-'p{t) near its 
maximum (rmax — 7.440, see above) and to the inverse of the angular frequency near the 
latter maximum ( ^ tt-fW — 7.968 ). 

^ <.^max 0.1255 ' 

In Fig. we display the temporal evolution of the recoil velocity, where we exhibit both 
its magnitude |v*Qjjj(t)| and its two separate components ^comi't) ^comi't)- We see that 
the maximum instantaneous recoil velocity is reached after the maximum of \J-'p\, and while 
|jFp| has already significantly decreased. After having reached its maximum value, |v*Qj^(t)| 
slightly decreases, in general with some oscillations, before settling down to its terminal value 
(see top panel in Fig. El). A useful visualization of the evolution of the recoil is provided by 
plotting the "hodograph" , i.e. a parametric plot of the instantaneous recoil velocity vector in 
the two-dimensional plane (v^^jj^, v^^^^): see Fig. IHl The behavior of the instantaneous recoil 
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velocity is easily interpreted in view of the analytical arguments presented above. Indeed, we 
have seen in Eq. that, during the plunge, the main contribution to the integrated recoil 
came from the non-adiabatic character of the evolution of |J^p(t)| near its maximum. This 
led us to estimate that the instantaneous recoil Vcom{t) was given by the truncated Gaussian 
integral, Eq. (|4ip. This integral can be expressed by a (complementary) error function 
erfc(2;) = (^-^^ 1^°^ with argument z = —\f\ where i = t — tmax, and 

where a and /3 are defined by Eq. (jl^ . If, for simplicity, one neglect e^ax ^ 1, we find that 

z ~ — ^ i^max Tmax j • Notc that z IS shlftcd (by a quantity of order unity because 

cUniaxTmax = Q — 0.9337 ) iu the complcx plane. This shift in the complex plane introduces 
some modifications to the usual behavior of the complementary error function in the real 
domain, which evolves monotonically from erfc(+cxD) = when z = +00, i.e. t = —00, to 
erfc(— 00) = 2 when z = —00, i.e. i = +00. These modifications are such that the modulus 
|erfc(2(t))| increases from the value when t = —00 to a maximum value of about 2.05 
when ~ +1.06, before decreasing toward its final value of 2 when ^ 1. Note also 
that |erfc(z(t))| already reaches the value ~ 1.97 when ^ +1. Therefore, most of the 
integrated effect of the maximum of |J?^p(t)| is acquired when t — tmax ~ +Tmax- However, 
in the case of the evolution depicted in Fig. |21 one can check that the time tmax + ^max 
corresponds to a radius r ^ 2.648, which is (slightly) below r ^ 3, i.e. after our chosen 
transition time to ring-down tmatch — 4153.5 (corresponding to r ~ 2.933). Therefore, the 
integrated effect up to tmatch of the non-adiabatic evolution of |^p(t)| near its maximum will 
be slightly smaller than the total integrated effect considered above. 

In addition, the transition from jFp"°^'^(t) to jFp°^(t) across tmatch introduces a new source 
of non-adiabaticity Fig. El shows that the ring-down behavior can introduce some oscilla- 
tions in |vcom|, and tends to decrease the value of |vcom| reached after passing the maximum 
of |^p|. However, one sees on the plot that the effect of ring-down is relatively small com- 
pared to the main contribution to |vcom| acquired by passing over the maximum of |jFp(t)|. 
The relative smallness of the ring-down contribution to |vcom| can also be checked ana- 
lytically. Indeed, this ring-down contribution is given by a sum of integrals of the form 
f^°° drC e~^°'~^^'^^'^ = p. . One can then relate this sum of integrals to the value of iJFpl 



The non-adiabatic character of the transition between plunge and ring-down shows up particularly in the 
fact that one passes from a quasi-monochromatic (chirping) form [J-p{t) + * •^p(0]piungc ~ a(i) e"^'*) to 
a non-monochromatic form containing both (decaying) positive and negative frequencies: e"'^e^*'^^. 
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at the moment of the matching. One then checks that, because the transition occurs while 
|jFp| has already significantly decreased from its maximum value, the ring-down integral will 
be significantly smaller than the value acquired by passing over the maximum of |jFp(t)|. 

In Figs. [7| and IHl we study the effect of demanding a smoother transition between plunge 
and ring-down, namely a one ( with two conjugate pairs of QNM's per multipole mo- 
ments) instead of the one ( one pair of QNM's) used in the Figs. |3] and El As we see, 
though the effect is not negligible (and can introduce some extra oscillations in v^^^^), it has 
a relatively minor effect on the final recoil velocity. More precisely, we find that (for rj = 0.2) 
vl~% QNM ^ 51.05 km/s, while vl~% qnm ^ 54.19 km/s. 

Let us now consider the impact on f*o^™'^^ of the higher-order PN corrections to jFp, i.e. 
the effect of the factor F in Eq. (jl8p . The definition of F, given by Eq. (j25|) . depends on the 
definition of F{vY'^. We have discussed above various ways of estimating F: one can use 
straightforward "Taylor approximants" , ( e.g. F'^pj^ = 1 + F2f^), or, instead, some "Fade" 
ones ( e.g. F[pjy = ^ r, \t] ) . Actually, the IFN-accurate Taylor approximant 

''polo ~^1 + C2V 

Fippf = I+F2 v"^ is not an acceptable approximation for the study of the recoil. Indeed, when 
Tj = 0.2, one finds that F^pj^f becomes negative for v > 0.421. As = u r is about 0.44 at 
the maximum of JFp, the most important domain of values for to estimate the recoil would 
correspond to such a physically incorrect negative value for F^pj^. Though the corresponding 
Fade approximant F^pj^ stays positive, it is also unphysical in that it takes values of order 
10 in the relevant range of values for v. Therefore, we shall only consider the higher-order 
FN approximants: F'^^p^^, F^p^, -^^lspat, and -F^ptv- TableHl we present (for 77 = 0.2), the 
values of |vcom|, v^om' ^com various approximants and at various stages of the inspiral 
or merger. We observe that changing the approximant for F has a substantial effect on Vcom- 
In particular, the final recoil varies between t'co™™'^^ — 50 km/s when using a 2FN accurate 
Fade approximant, and fco™"^' — 74 km/s with a 2FN accurate Taylor approximant (the 
quasi- Newtonian estimate being f*™™'^^ ^ 51 km/s). In agreement with the approximate 
analytical estimate derived above, one can check that Wco™°^^ varies in direct proportion to 
the value taken by the FN factor F at t = tmax, i-e. when |J?^p(t)| reaches its maximum 
value. Indeed, neglecting the effect of ip in Eq. [^(^max) ~ 1.007], the variation of 



For simplicity, in view of the closeness of ijj{r,pip) to 1, we shall not contemplate other definitions of F 
based, e.g., on expanding both ip and F in PN series before, eventually re-summing the PN expansion of 
F. 
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'^co™'^'^^ with F is well describable by the variation of F{vma.^ with f max — 0.439. Indeed, we 
have the following values for F, namely, -F/v(fmax) = 1, F]-^p^{vraa,x) — 1.33, Fj^ 



2PAfV^maxj — 

1.35, Ff^pj^iv^i^^) ^ 1.27, and -F2^jv('^max) — 0.93 which are well correlated with the results 
listed in Table H] The significant difference between -F2'pAf(^max) and -^2^7v(^'max) illustrates 
again the poor convergence of the successive PN contributions. As argued earlier, one 
would expect [by analogy with convergence of the energy-flux function FEiv)] that the 2.5 
PN-accurate Pade approximant would yield a better answer. In absence of information 
concerning the 2.5 PN level, we shall use -F2PAr hesi answer [in view of the behavior 

of -^Eap^ but keep in mind the probability of a significant error bar around it. 

In Table IHl we study the influence of another parameter in our methodology: the precise 
choice of the transition radius between plunge and ring-down. We consider the standard 
Schwarzschild light ring, namely r = 3, as our default value. In Table m we explore the 
effect on ^'0™"^' of changing the matching radius rmatch by ±20% of its default value. As 
this Table shows, the value of fco™™'^^ is mildly sensitive to the precise choice of transition 
radius. 

Up to now, we have only focused on a specific symmetric mass ratio (77 = 0.2) which, 
in view of the analytical estimate, given by Eqs. (jl^ and (jlHI), is expected to yield the 
maximum possible recoil. In Table IIIH we consider several different values of 77, namely 
7] = 0.05, 0.1, 0.2, and 0.24 and compute the corresponding sca/ec? terminal recoil Vcom = ^f^^ 
where /(r/) = t]'^ ^/T^^Tr]. As expected, after scaling out the function /(//), the recoil 
depends only weakly on rj. We can analytically approximate the ?7-dependence of Vcom by a 
second-order polynomial P{ri) = a + hrj + erf . Normalizing a, h and c so that P(0.2) = 1, 
we find a reasonable fit^^ for 

P(r/) = 1.0912 - 1.04 7/ + 2.92 r/^ (59) 

Finally, putting together the various informations we have obtained above, we can sum- 
marize our "best bet" estimate for the final recoil associated with the coalescence of binary 
black holes of symmetric mass ratio 77 as 



^com ^73.5^/(r/)km/s, (60) 



Actually, to get a good fit to the data in Table UTTl one needs a third-order polynomial, namely: P-iitf) 
1.112 - 1.78 + 10.3 if - 21 if. 
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where /(r/) = f^ff^ (1.0912 - 1.04 r/ + 2.92 r/^). The fiducial value 1.35 used above for 
scaling Fmax is the prediction made by the 2PN Taylor approximant around rmax — 3.50, i.e. 
at the moment where the modulus of the linear momentum fiux is maximum. Note that the 
proportionality of the final result to -Fmax is only approximate because the presence of the 
correction factor F(t) changes not only the height of the maximum of |jFp| but also affects 
the shape of |J^p|(t) and thereby the quality factor Q etc. The above estimate is plotted as 
a function of g = — in Fig. El 



VI. DISCUSSION 



The main fruit of the present study is the fact that we have delineated, often by means 
of analytical arguments, the relative importance of several different physical effects in de- 
termining the magnitude of the final recoil velocity fcom- We have emphasized that the 
value of f c,5m^ is essentially determined by a brief period during the orbital evolution when 
the integrand of the oscillatory integral ()34|) yielding f com varies in a non adiabatic manner: 
a/a ~ 0. We have found that this non-adiabatic evolution is confined to a small neighbor- 
hood of the moment, during the plunge, where the modulus of the linear momentum fiux 
|jFp|(t) [i.e. the amplitude a(t) in the integral (jS3)] reaches a maximum. The good news is 
that it seems that this maximum takes place during the quasi-circular "plunge phase" , i.e. 
during a phase where the radial kinetic energy is significantly smaller than the azimuthal 
kinetic energy. Indeed, the ratio TZ = g^^p^/ g'^'^p^^p between "radial" and "azimuthal" ki- 
netic energies is found to take the value T^max — 0.135 at rmax- [Note again, that even at 
the light ring, r ^ 3 this ratio remains small, namely, TZ\r ^ 0.281]. 

This "burst" of linear momentum fiux also occurs slightly before the merger phase which 



we view as taking place when the^ 
about 3. As was argued in Ref. 



adimensionalized) radial distance r gets smaller than 
I the quasi-circular "plunge" phase (3 < r < 6) is a 
priori amenable to analytical description within the EOB approach. And we have indeed 
checked that various different ways of completing the EOB approach by a suitable matching 
to a subsequent (r < 3) ring-down description of the merging of the two black holes did 
not affect much the recoil velocity acquired after passing over the maximum of |jFp|(t). 
We have also verified that various other physical ingredients of the model (such as: the 
representation of the damping force during the plunge, the choice of matching point, the 
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number of quasi-normal modes, • ■ ■ ) had a rather mild effect on the final recoil. 

However, the bad news is that when |J^p|(t) reaches its maximum value, the azimuthal 
kinetic energy contribution p^/r^ in the Hamiltonian equations, Eqs. (jl7p . as illustrated 
in Fig. ^ is of order unity ( p^/r^ax ~ 1)' i-^- comparable to the constant term (= 1), 



which plays the role of the "squared rest mass" term in the Hamiltonian {H ~ y + m^). 
This situation contrasts with the one near the LSO, where one has p'^/rlso ~ 0-25 which 
is significantly smaller than unity. In other words, the orbital motion near the LSO is still 
"non relativistic" (by a thin margin), while the recoil is generated when the orbital motion 
becomes mildly relativistic (in the sense ~ m^). What further complicates the matter is 
that the orbital motion does not follow the usually considered sequence of circular orbits, so 
that we cannot use the standardly assumed relativistic version of Kepler's third law relating 
the angular velocity to the radius. In the body of the paper, we have used the "quasi- 
Newtonian" expression for the linear momentum flux (in EOB coordinates) as a guideline 
to select a "best bet" modeling of J-'p(t) during the plunge. We have already seen in Table HI 
that a very significant source of uncertainty in the magnitude of J^p(t) concerns the inclusion 
of post-Newtonian corrections in it. Depending on whether one uses the straightforward 
"Taylor-expanded" 2PN correction ^j, or one of its Pade-resummed versions, one gets a 
multiplicative factor varying between 0.92 and 1.35 in v^^^. However, this uncertainty is only 
a lower limit to the total uncertainty currently attached to the description of the relativistic 
effects during the plunge. We can see hints of a larger uncertainty by comparing our EOB- 
based treatment (which did not assume the validity of the standard Kepler law during the 
plunge) to a treatment similar to the one advocated in Ref. (which did implicitly assume 
the continued validity of Kepler's law). To explore this issue, and also to understand the 
relation between our estimate and the significantly larger one obtained in Ref. P], we have 
estimated the recoil following from using the functional form (jS)) for jFp, instead of our 
2PN-corrected, quasi-Newtonian expression (|T8|l . As we discussed above, these two different 
functional forms can be matched above the LSO [modulo the suitable definition of the 2PN- 
correcting factor F in Eq. ()18|1 . see Eq. ()25|l ]. by using the relativistic Kepler law ()2H) . On 
the other hand, Eq. (PT|) gets strongly violated below the LSO, and as a consequence the two 
different functional forms, Eqs. Q and ()18p. a priori lead to quite different time evolutions 
for JFp(t). In keeping with our general philosophy, it is useful to understand analytically the 
difference between the two basic corresponding "quasi-Newtonian" prescriptions [obtained 
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by neglecting the 2PN factor F in Eq. Q and the corresponding 2PN factor F in Eq. p8|) ]. 
In other words, let us discuss the effect of replacing our basic "fiducial" quasi-Newtonian 
momentum flux !Fr,u) = f'^'jJ^ by T^^ = v^j^ = (cuY^^'^. Let us call K = u'^r^ the quantity 
which is (approximately) constant when Kepler's law is satisfied. It is then easy to see 
that the ratio between the two prescriptions reads: T^^jTr^uj = K~^/^. We have shown 
above how to write an approximate evolution equation for the angular frequency during the 
plunge, namely: uj oc A{r)/r'^. This entails that K varies during the plunge approximately 
as oc A^{r)/r. It can be analytically checked that or better K = ipK with ip defined by 
Eq. (j22I) (augmented by the needed pi terms), is equal to 1, and has a horizontal derivative, 
at the LSO, and decreases monotonically when r < tlso, fo reach zero when r tends to 
the "horizon" [A{r horizon) = 0]. This behavior is illustrated in Fig. 01 As a consequence K 
gets (significantly) smaller than 1 below the LSO, so that Tv^I^t.uj = K'^l"^ is significantly 
larger than 1 during the plunge. More precisely, we can, as above, study the evolution of 
Tv^{t) by writing that it approximately varies like = [uj{r)Y^^'^ oc [A(r)/r^]-'^^/'^. It can 
be seen that this function of r has a maximum at r^g^^ ~ 2.879, and that the value of the 
function at its maximum is more than twice higher than the maximum that J^r,w = r^uj{ry 
had at rmax ~ 3.443. This increase in the maximum value of the linear momentum fiux 
is further compounded by significant changes in the shape of the maximum (notably the 
value of the crucial quality factor Q = cUmaxTmax introduced above). We therefore see that 
using a momentum fiux given by Eq. ^ instead of Eq. (fT^ will more than double the value 
of f^oj^'[see precise numbers below]. Note that the large change in the predicted value for 
v^^^^ that we just discussed concerns only the "leading-order quasi-Newtonian" expression 
that one chooses to employ during the plunge. It has to be further compounded by the 
uncertainty linked to the resummation of the 2PN correction factor F ot F (which brings a 
multiplicative uncertainty factor of order 1.5). 

Actually, the fact that the location of the maximum of the momentum fiux (r^^^^^ ^ 2.879) 
is slightly below the light ring in the case of Eq. (jSj) brings a further complication. Indeed, in 
that case we cannot trust our simple "Gaussian integral estimate" Eq. ()44j) which assumed 
that one was integrating over the maximum of |jFp|(t). The matching to the subsequent 
ring-down behavior could a priori have a significant impact. To resolve this issue we have 
run numerical simulations similar to what we use in the body of the paper (with EOB 
dynamics as defined above, and two Quasi-Normal-Modes matching at r = 3) except that 
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the right-hand side of Eq. (fTSj) was replaced by Eq. Q. Our results are displayed in Fig. [TUl 
and listed in the last two rows of Table HI In agreement with the simple analytic arguments 
above, we do find final recoil velocities that are more than twice larger when using Eq. (|18p 
(with a corresponding factor F). For instance, from Fig. ^[ we find that the terminal 
recoil is ~ 172km/s which is more than double the value that we read in the Table U for a 
2PN accurate F. Furthermore, using Eq. Q as it stands we observe, for the optimal case 
7] = 0.2 and while terminating the EOB evolution at r ~ 2 without doing QNMs matching, 
"^com^ — 243km/s. On the other hand, if we do not include the 2PN correction factor F 
in Eq. (jH)) we obtain, again for t] = 0.2, a final recoil ~ 135km/s at r ^ 2. These results 
are roughly consistent with the results of Ref. P], and confirm our diagnostics that the 
main physical origin of the integrated recoil is the rather well-localized "burst" in linear 
momentum flux occurring during the plunge. 

One might view the large difference between the two models of momentum flux, Eq. p8|) 
versus Eq. 0, in several different ways. One way would be to say that, as we have seen, it 
is not justified to continue assuming (as is implicitly done in Eq. ^) the validity of Kepler's 
law during the plunge, and therefore that the corresponding prediction for Vcom is definitely 
too large^^. On the other hand, we have emphasized above that near the crucial maximum 
of radiation of linear momentum, the orbital motion becomes mildly relativistic (p^ ~ m^, 
i.e. {v/cY ~ 0.5) and, in addition, more complicated than the quasi-circular and quasi- 
Keplerian cases studied so far in detail in analytical gravitational wave research. Therefore, 
one might also say that the large difference between the two proposed extrapolations for 
the momentum flux just reflects our ignorance of what is the correct momentum flux in 
such a relativistic situation. We do tend to think that our "best bet" estimate, Eq. (lUUj). 
is probably closer to the truth, but we cannot provide any proof of this belief, nor can we 
presently define an "error bar" around our preferred estimate fj60p . One way to estimate an 
"error bar" around (j60|) would be to study the effect of using a 3PN-accurate EOB dynamics, 
and/or to include all the contributions proportional to R and R (which were neglected here) 
in the momentum flux. On the other hand, we consider it likely that the results quoted 
above (and listed in the last two rows of Table coming from the Kepler-law based Eq. 

Because, as we have just seen, the known violation of Kepler's third law, i.e. K < 1, is the root of the 
difference between the two estimates. 
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furnish an upper bound on the correct recoiL 

We conclude that none of the current analytical-based estimates of the total recoil are 
reliable. From Table El we see that, depending on the analytical representation used for 
the linear momentum flux during the plunge, we get a final recoil velocity which varies in 
the range 49 — 172 km/s. This "uncertainty range" in theoretical predictions is illustrated 
in FigEl We view this current theoretical uncertainty as a strong motivation for devoting 
future work to the specific issue of selecting a reliable^^ analytical model of linear momentum 
flux during the plunge. There are several avenues one could use towards this aim. One might 
generalize the work of Ref. by keeping the terms proportional to or p^, and by making a 
minimal use of Kepler's third law^^. Another avenue is to compare analytical and numerical 
results in the case of a test-particle, 77 <^ 1, plunging into a Schwarzschild black hole. One 
might also try to compare analytic and numerical results in the case of full 3-d simulations of 
coalescing black holes. However, one should keep in mind that the linear momentum flux 
is a sub-dominant effect in the gravitational wave emission which can easily get drowned in 
the "noise" associated, for instance, with the presence of residual incoming radiation in the 
initial data. Let us note in this respect, that the situation is a priori much better in the 
(more urgent) problem of the modeling of gravitational radiation (linked to the dominant 
energy flux) from coalescing black holes. In this case most of the signal-to-noise ratio is 
linked to the train of waves emitted in the last few orbits before crossing the LSO. The EOB 
formalism was conceived for dealing with this problem, and the current study should not be 
interpreted as casting any serious doubt on previous studies j^, [ill which relied mainly on 
rather robust analytic features of the EOB approach. 



By "reliable" we mean here "accurate within 50 %" or so. The current uncertainties in analytical estimates 
span a factor ^ 5. 

For instance, one might use Kepler's law only at the level of the derivation of the multipole moments, and 
not use it anymore when time-differentiating them along plunging orbits. One might also use resummation 
techniques, notably by using the resummed quasi-Schwarzschild coordinates that enter the EOB formalism. 
After the submittal of this paper, remarkable progress in numerical relativity has allowed Baker et al. 
[45! to report on the first accurate numerical calculation of the recoil velocity of two non-spinning black 
holes with 77 = 0.24. Their result (wrocoii = 105 km/s) corresponds, if we use the scaling function f{i]), 
defined after Eq. (|60|l . to a maximum recoil of 161.5km/s. This value falls within the range of analytical 
estimates summarized in Table More precisely, the last row of Table ^ shows that the momentum flux 
Eq. Q predicts a recoil of 110 km/sec when 77 = 0.24. We note also that the time evolution, during 
the coalescence, of the magnitude of the instantaneous recoil velocity showed in Fig. 1 of Q is in good 
qualitative agreement with our analytical predictions (see Figs. 151 151or ll0() . 
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If our 'best bet' estimate, Eq. ()6Up. is confirmed, it might have significant astrophysical 
consequences. For instance, a low value for the terminal recoil will have implications for the 
formation of massive black holes at high redshifts j^l, which may influence various event 
rates for binaries involving supermassive black holes that LISA may observe. Finally, a recent 
surge in astrophysical investigations that probe consequences of the recoil indicates that it 
is desirable to know the dependencies of recoil on the orbital configurations of coalescing 
black hole binaries ^|. Therefore, it will be also important, in the near future, to extend 
our approach to the study of the recoil associated with the coalescence of spinning black 
holes in inspiralling eccentric orbits. 
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TABLE I: Values of I 

Vcomli v^Qjjj and Vcoiri) recoil velocity and its x and y components in kms 
for 7] = 0.2 at various stages of coalescence for the momentum flux Eq. (|18|) with different F's. 
The last two separated rows (corresponding to rj = 0.2 and 7] = 0.24) show the recoil estimates 
obtained by using, instead of Eq. p8|) . the (Kepler-law-based) momentum flux Eq. 





rLSO - 6.00 


y-match =i 3.00 


t — > -|-oo 


1.5 PN Taylor F 


22.09, 18.05,12.73 


82.84,51.73,-82.68 


72.41,12.18,-71.37 


2 PN Taylor F 


22.34,18.25,1.29 


84.07, 5.23,-83.91 


73.52,12.31,-72.48 


F = 1 


18.50,15.41,10.24 


59.83,4.29,-59.67 


51.05,10.29,-50.00 


1.5 PN Pade F 


20.69,16.91,11.92 


79.97,5.11,-79.81 


70.19, 11.68,-69.21 


2 PN Pade F 


16.20,13.42,9.08 


56.68, 3.96,-56.54 


49.02,9.15,-48.16 


Eq. (H, 7? = 0.2 


23.39, 19.34,13.16 


179.88, 78.54, -161.82 


171.55, 136.64,-103.73 


Eq. (H, r? = 0.24 


15.47, -5.03,-14.62 


118.67,-85.74,82.05 


109.64,-103.45,36.34 



TABLE IL Effect of changing the transition point where merger phase goes to QNM ringing on 
terminal |vcom| for t] = 0.2 and F = 1, when changing the canonical value, r^^^ch — 3, by ±20%. 





1 Vcom 1 


* com 


Vcom 


'"match — 3.6 


44.98 


-8.05 


-44.26 


'■match ^ 3.00 


51.05 


10.29 


-50.00 


''match — 2.4 


58.21 


12.60 


-56.83 



TABLE IIL Dependence of scaled terminal recoil velocity ^2 on r] for rmatch — 3.00 and 

F =1. 





In kms ^ 


7] = 


0.24 


2817.63 


7] = 


0.2 


2853.81 


7J = 


0.1 


2901.29 


7] = 


0.05 


2987.44 
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FIG. 1: Plots of 'azimuthal' and 'radial' kinetic energies as functions of coordinate time t and 
radial separation r during the late inspiral and plunge for a i] = 0.2 binary. The initial orbital 
separation, when t = 0, was r = 15, and the curves were terminated at t ~ 4153, corresponding 
to r = 3. The plots show that dominates pr{t)'^ / B{r) during the entire EOB evolution, 

including the plunge (down to r = 3). 
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FIG. 2: Magnitude of the linear momentum flux (in the quasi-Newtonian approximation) as a 
function of r (left panel) and t (right panel) during the late inspiral and plunge for a r/ = 0.2 
binary, whose initial orbital separation, when t = 0, was r = 15. We terminated the plunge 
arbitrarily around r ~ 2.65. 





42 



11=0.2 




FIG. 3: Plots of lo, uP' and , / iii terms of r resulting from the EOB evolution, given 

by Eqs. (|17() . for r] = 0.2 binary. The panels clearly demonstrate that the orbital frequency evolves 
differently during the late inspiral and the subsequent plunge. Note in particular the strong decrease 
of the "Kepler combination" K = r'^ up' during the plunge. 
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FIG. 4: Plots for the magnitude, x and y components of the linear momentum flux (quasi- 
Newtonian approximation) versus t during the coalescence of a binary with r] = 0.2. For this 
figure, the ring-down phase is described by 2 QNMs and orbital separation of the binary, when the 
numerical evolution began, was r = 15. 
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FIG. 5: Temporal plots for the magnitude, x and y components of recoil velocity for the binary 
configuration described in Fig. = 0.2, quasi-Newtonian flux, 2 QNM's). 
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FIG. 6: Parametric plot of ^com versus Vcom during the coalescence of a = 0.2 binary. We 
employ the quasi-Newtonian approximation to the linear momentum flux followed by a 2 QNMs 
description of the ring-down phase. The symbols x and * respectively indicate the positions where 
the linear momentum flux reached its maximum value and where matching to the ring-down phase 
was done. The initial orbital separation of the binary, when f = 0, was r = 15. 
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FIG. 7: Plots for the magnitude, x and y components of the quasi-Newtonian hnear momentum 
flux versus t during the coalescence of a binary with rj = 0.2. In this case, the ring-down phase is 
described by 4 QNMs and the initial orbital separation of the binary, when t = 0, was r = 15. 
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FIG. 8: Temporal plots for the magnitude, x and y components of recoil velocity for the binary 
configuration described in Fig. [TK?? = 0.2, quasi-Newtonian flux, 4 QNM's). 
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FIG. 9: Our 'best-bet' estimate (solid line) for the terminal recoil as a function of the mass 
ratio 9 = ^ based on Eq. (|6U() with F = 1.35, corresponding to using the 2PN, Taylor F in the 
momentum flux expression Eq. ()18() . The "theoretical uncertainty" around this best-bet estimate 
is illustrated by plotting the results of using: (i) a quasi-Newtonian flux (Eq. (|18j) with F = 1) 
[dot-dashed lower curve], or (ii) the Kepler-law-based 2PN flux Eq. ^ [dotted upper curve]. The 
estimates taken from Tables HI and IIIl| are denoted by + symbols. 
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FIG. 10: Temporal plots, along the EOB evolution, for the magnitude of the "Blanchet-Qusailah- 
Will" linear momentum flux, defined by Eq. (j^J, the orbital separation r and the associated recoil 
I f com I in km/s for rj = 0.2 binary. We terminate the plunge around r ~ 3 and perform two QNMs 
matching for the ring-down phase. [In Ref. the evolution was formally continued down to r ~ 2]. 
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